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ABSTRACT 

We investigate the effect of ionising radiation from Massive Young Stellar Ob- 
jects impinging on their emerging spectral energy distribution. By means of detailed 
radiative transfer calculations including both the gaseous and dust phase of their sur- 
rounding protoplanetary discs we highlight the importance of modelling both phases 
simultaneously when interpreting observations from such objects. In particular we find 
that models that only include dust may lead to incorrect conclusions about the inner 
disc evolution. Furthermore the omission of gas from models overproduces far-infrared 
and sub-millimitcr fluxes with the result that derived dust masses may be underes- 
timated by a factor of two in some cases. Finally free-free emission from the ionised 
component of gaseous discs causes the slope of the dust emission in the sub-mm and 
mm regime to appear flatter, resulting in incorrectly modelling the dust properties, 
with consequences on the derived disc masses, power law index of the surface density 
profile and other disc properties. 
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1 INTRODUCTION 

Over the past two decades, the number of observations of 
protoplanetary discs has rapidly accelerated, and with it our 
understanding of the properties of these objects. With ad- 
vances in resolution and precision, instruments such as the 
Sub-Millimeter Array (SMA) and Spitzer have allowed for 
the detection of many hundreds of protoplanetary discs in 
various environments. In particular, large surveys performed 
by Spitzer have mapped ~ 90% of all the star-forming re- 
gions within 500pc of the Sun and have obtained spectra for 
over 2000 YSOs contained therein. 

Observations of those star-forming regions nearest to 
us have, to date, dominated the literature as they con- 
tain predominantly low mass, comparatively isolated stars 
which allow the structure of discs and their evolution to 
be determined more easily than their higher mass counter- 
parts. These include regions such as Taurus, Ophiucus and 
Chamaeleon I. More recently, studies have looked slightly 
further afield to the Orion nebula cluster where, in addition 
to studies in the infrared, the HST has beautifully captured 
images of protoplanetary discs in the optical. 
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To date, nearly all identified protoplanetary discs sur- 
round low or intermediate mass stars, with only a few de- 
tected around stars with masses > 8Mq. Despite many 
searches, no convincing evidence for the existence of ac- 
cretion discs around protostars with masses > 20Mq (O- 
type stars) has yet been found (e.g., Cesaroni et al. 2007). 
Moreover, only a very limited number of potentially disc- 
harbouring sources with masses > IOMq (early B-type 
stars) are known. One of these is W33A with a mass of 
10 to 15M0 (Davies et al., 2010). A non-spherical, compact 
component was detected with MIDI at the VLTI (de Wit, 
2010). Another example is IRAS 13481-6124 with a mass 
of about 20Mq. Kraus et al. (2010) reported from interfer- 
ometric measurements in the near-infrared, obtained with 
AMBER at the VLTI, a compact, elongated, hot inner com- 
ponent. A final example is the Silhouette Disk around a 
young star in M17 (Chini et al. 2004). The mass of the cen- 
tral object is uncertain, however, and could be as low as 3Mq 
(Nielbock et al., 2008). No circumstellar disc surrounding a 
star larger than ~ 20AfQ has so far been reported (Preibisch 
et al 2011). 

The existence of circumstellar discs around massive 
young stellar objects (MYSOs) is still sometimes questioned 
in the literature, as it remains unclear by what process high- 
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mass stars form. Frequency dependent radiative feedback 
simulations by Kuiper et al. (2010) indicate that it is possi- 
ble for high mass stars to form by the same mechanism as 
lower mass stars. Hovsrever, it is not clear how reliable these 
results are as the numerical resolution of the simulations is 
low (1.27AU). An alternative to this process is discussed by 
Zinnecker and Yorke (2007) whereby several low-mass stars 
merge together to form a high-mass star. 

Even assuming the existence of discs around MYSOs, 
the detection of these objects is a significant challenge: mas- 
sive stars generally form at large distances, they are nor- 
mally deeply embedded in the cloud in which they form and 
they evolve very quickly making them relatively rare. Ad- 
ditionally, the discs themselves have a very short dispersal 
timescale and thus MYSO-disc systems are likely very rare. 

However, very recent studies by Preibisch et al. (2011) 
and Grellmann et al. (2011) suggest the discovery of the po- 
tential first discs around high-mass stars. Both focus on a 
single disc-star system, the first in the Carina Nebula and 
the second just north of the Cone Nebula. Both studies fit 
the model SEDs from Robitaille et al. (2006) to their obser- 
vations, and derive approximate stellar masses of ~ 12Mq 
and ~ 11.6AfQ respectively. 

This technique is valuable as it allows information about 
the structure, composition and evolution of discs to be dis- 
cerned. Heretofore, models used for this purpose have only 
needed to simulate stars at lower masses and temperatures, 
generally of A/F type or later, since observations were gener- 
ally restricted to stars of these types. At these temperatures 
(^ 7000K), the stars have very little or no ionizing compo- 
nent to their energy spectrum (assuming a blackbody energy 
distribution). As a result, the models for discs around low 
mass stars need only to consider radiative transfer via the 
dust particles in the discs, as the gas opacities are orders 
of magnitude lower than the dust opacities at non ionising 
frequencies. 

However, the new discovery of discs around MYSOs cre- 
ates a need for more accurate models. Future observations 
will hopefully allow for the study of more discs around even 
more massive YSOs (> 20Mq), where a significant fraction 
of the stellar continuum is emitted longward of the Lyman 
edge. For this reason, it may be necessary to consider the 
effect of interactions between the stellar ionizing photons 
and the gas component in the disc on the spectrum that 
the dust is exposed to, and thus on the overall SED of the 
disc-star system. With a typical disc gas-to-dust mass ra- 
tio of ~100, it is likely that the ionizing component of the 
star's radiation field will interact with gas particles before 
being absorbed by any dust particles, and thus that the dust 
will be exposed to an overall softer radiation field than in 
the "dust-only" case. This could potentially result in less 
emission than expected from the dust in the mid-IR region, 
and thus lead to an underestimate of the mass contained 
within the disc. This is particularly relevant if there is any 
significant difference between prediction of models with and 
without gas at wavelengths that are frequently used for ob- 
servations; we will consider differences in the models at the 
Spitzer Infrared Array Camera (IRAC) wavelengths (3.6, 
4.5, 5.8 and Sfim) and the Spitzer Multiband Imaging Pho- 
tometer (MIPS) wavelengths (24, 70 and 160/im), in addi- 
tion to those wavelengths used to observe the circumstellar 



disc surrounding the Carina Nebula (see Table 1 (reference 
1), Preibisch et al. 2011). 

In section 2 of this paper we describe our numerical 
approach and modelling strategy, in Section 3 we present our 
results, which are then discussed in more detail in Section 
4. Finally our main results and conclusions are summarised 
in Section 5. 



2 NUMERICAL APPROACH 

We use the 3D photoionisation code MOCASSIN (Ercolano 
et al., 2003, 2005, 2008a) to produce the SEDs of disc- 
star systems for stars at temperatures of 5800K, lOOOOK, 
20000K, 30000K, 40000K and 50000K. The code uses a 
Monte Carlo approach to the transfer of radiation, allow- 
ing the treatment of both the direct stellar radiation as well 
as diffuse fields and the transfer through dust and gas phase. 
The code includes all the dominant microphysical processes 
that infiuence the gas ionisation balance and the thermal 
balance of dust and gas, including processes that couple 
the gas and dust phases. In the ionised region the dom- 
inant heating process for typical gas abundances is pho- 
toionisation of hydrogen, which is balanced by cooling by 
coUisionally excited line emission (dominant), recombina- 
tion line emission, free-bound and free-free emission. The 
atomic database includes opacity data from Verner et al. 
(1993) and Verner & Yakovlev (1995), energy levels, colli- 
sion strengths and transition probabilities from Version 7 of 
the CHIANTI database (Landi et al. 2006, and references 
therein) and the improved hydrogen and helium free-bound 
continuous emission data of Ercolano & Storey (2006) . The 
code was originally developed for the detailed spectroscopic 
modelling of ionised gaseous nebulae (e.g. Ercolano et al 
2004, 2007), but is regularly applied to a wide range of astro- 
physical environments, including protoplanetary discs (e.g. 
Ercolano et al 2008b, 2009, Owen et al 2010, Schisano et al 
2010). Arbitrary ionising spectra can be used as well as mul- 
tiple ionisation sources whose ionised volumes may or may 
not overlap, with the overlap region being self-consistently 
treated by the code. Arbitrary dust abundances, composi- 
tions and size distributions can be used, with independent 
grain temperatures calculated for individual grain sizes, al- 
lowing to self-consistently calculate the sublimation radii of 
grains of different sizes. 

We use blackbody input spectra for the stellar radia- 
tion at the effective temperatures specified above and for 
each temperature we consider four different optical depths 
(where we define the optical depth at the midplane of the 
disc Tmp) varying from the optically thin model Tmp = 0.1 
to the optically thick models Tmp = 1.0, 10.0 and 100.0. We 
describe the dust density distribution similarly to Pascucci 
et al. (2004) by: 
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where a = 1, and /3 — 1.125, leading to a surface density 
profile of E(r) cx ^--o i^s^ ^]^g Q^jjgj- (;jjgg parameters are de- 
scribed in Table 1 of Pascucci et al (2004) and include a disc 
inner radius of 1 AU and a disc height of 125 AU. 

We use the dust density distribution to compute the 
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density distribution of the gas, setting the gas-to-dust mass 
ratio equal to 100. We discuss the results along two different 
lines of sight through the disc: one nearly edge-on {i —77.5°) 
and the other nearly face-on {i =12.5°), where i is the incli- 
nation of the disc. 

We run two sets of models, one with only dust included 
in the grids and the other where both dust and gas are 
present. All models where run for solar metallicities. We con- 
sider the percentage flux differences between the "dust-gas" 
cases and the "dust-only" cases. In particular, we consider 
the flux differences at the observational wavelengths used 
in Preibisch et al 2011 in addition to those wavelengths 
used by Spitzer (IRAC & MIPS) (see Table 1 for a sum- 
mary). Each model is defined by the optical depth, stel- 
lar temperature and line of sight and is labelled using the 
following convention "t[optical depth]-T[stellar temperature 
in 1000K]-i[MP/POL]", where MP indicates the near mid- 
plane line of sight (77.5°) and POL indicates the near pole 
line of sight (12.5°). For example, the model with parame- 
ters r„p = 10.0, Teff = 20000K and i = 12.5° is labelled as 
"tl0.0-T20-iPOL" . The percentage difference is given as a 
percentage of the "dust-only" flux. 



3 RESULTS 

Table 1 shows the percentual difference between the dust- 
gas and dust-only model fluxes at specific wavelengths for all 
models with Teff > 20kK. We also include the 24 SED plots 
for each model in Figure [T] (for r,„p — 10, as a case study) 
and in Appendix A (for Tmp ~ 0.1, 1, 100). We will discuss 
the key differences between the two sets of models according 
to the appearance of the SEDs at specific wavelength ranges, 
specifically at wavelengths < O.l^m, in the range ~0.1-6/^m, 
in the range ~6-300/xm and at wavelengths> 800/xm. As ex- 
pected, we see no differences between the dust-gas case and 
the dust-only case at the two lowest effective stellar tem- 
peratures considered {T^s = 5800/C and T^s = 10000/C), 
with differences first appearing at Tcfi — 20000/(', and be- 
coming larger for larger stellar temperatures and with larger 
optical depth. This is easy to understand since only mod- 
els with Tefi J? 20000-/? produce enough photons at energies 
greater than the hydrogen ionisation potential. While these 
low Toff models served their purpose of benchmarking the 
gas and dust transfer against the dust-only transfer, they 
are not useful in the context of the scientific discussion to 
follow and will not be further considered. 

We now consider the potential effects of ionizing pho- 
tons on both the dust and the gas in the disc with reference 
to the differences between the dust-gas and dust-only SEDs. 
We assume in our models a gas-to-dust mass ratio of 100, 
which results in a much higher number density of gas parti- 
cles than dust particles, thus a significant percentage of the 
ionizing photons are absorbed by the gas particles before 
they have any interactions with the dust particles. If this is 
the case then more ionizing stellar photons are absorbed by 
the disc in the dust-gas model compared to the dust-only 
model. This effect can be indirectly seen in Figure [T] as a 
reduced flux in the edge-on dust-gas case for wavelengths 
< 0.09/im, which is due to the extinction of the light from 
the star along the line of sight to the observer. At all op- 
tical depths, this discrepancy can only be seen at temper- 



atures J5 20000-/?. Since the pole-on line-of sight does not 
go through the disc, the effect is not seen for that viewing 
angle. 

In both dust-only and dust-gas cases we see thermal 
emission in the infrared by the dust. Since, as was mentioned 
above, a significant percentage of the stellar ionizing photons 
are absorbed by the gas particles rather than dust particles, 
the spectrum that the dust is exposed to is, overall, softer 
in the dust-gas models than in the dust-only models. Part 
of the stellar flux absorbed by the gas will be reprocessed to 
less energetic (ionising or non-ionising) continuum or lines. 
The total flux as a result of dust emission would therefore be 
reduced in the dust-gas case relative to the dust-only case, 
and this accounts for the differences we see in the SEDs in 
the wavelength range ~10 - 300/-im. 

Given that the gas is absorbing energetic photons, we 
expect to see some evidence of gas emission, which we do 
in the ~0.1 - lO/im wavelength range. The shape of the 
curves are different in this range, with the dust-only case 
exhibiting a smooth curve whilst the dust-gas SEDs has the 
typical "step-like" appearance due to free-bound processes 
that dominate the gas emission coefficient at these wave- 
lengths. As one may expect, this feature becomes more pro- 
nounced with both increasing optical depth and temperature 
see Figures [U EH El andlM]). We note that Robitaille et 
al. (2006) found that the IRAC and MIPS colors of discs 
around massive stars for dust-only models would be very 
red - the hottest dust in the disc is set by the sublima- 
tion temperature, but the stellar temperature can increase, 
causing a gap in wavelength space (this can be seen e.g. in 
Figure [T] where the bottom panels have large gaps between 
the photospheric and disc emission for the dust-only mod- 
els). Taking into account the emission from the gas in the 
near- and mid-IR mitigates this effect, so that the very red 
colors predicted for discs by Robitaille et al. may not be seen 
in real MYSOs. 

Finally, since a fraction of the stellar blackbody spec- 
trum is ionising, we expect some of the gas particles to be 
ionised and thus for there to be free electrons present in 
the disc. The flux seen at long mm wavelengths (lOOO^m < 
A < lOOOO/im) in the dust-gas case is indeed due to free-free 
emission from the ionised gas component of the disc. 

As the effective stellar temperature increases, so does 
the ionizing fraction of the stellar continuum. It is then log- 
ical to expect that all of the above effects would therefore 
be enhanced with increasing Tea, which is why the flux dif- 
ferences between the dust-gas and dust-only cases become 
more pronounced the higher the temperature of the central 
star. 

To summarize these effects, we have shown in Figure [2] 
the ratio of the dust-gas SEDs to the dust-only SEDs. 



4 DISCUSSION 

Multiwavelength photometry may be used to construct 
SEDs of observed disc-star systems. By using fitting tools 
such as that detailed in Robitaille et al. (2007), a set (or 
sometimes multiple sets) of parameters may be identified 
that reproduce the observations and thus provide an in- 
sight into the physical properties of the disc-star system. 
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Figure 1. SED plots for mid-plane optical depth Tmp = 10. From left to right, the SED plots are for the near-pole line of sight (12.5°), 
the near mid-plane line of sight (77.5°), and averaged over viewing angles respectively. From top to bottom, the SEDs are for stellar 
effective temperatures of 5800K to 50,000K. Each panel shows the SED for the dust-only case (black), and the dust-gas case (red to 
blue). SED plots for the Tmp = 0.1, 1, and 100 cases are included in Appendix [XI 
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Table 1. Difference between dust-gas and dust-only flux as a percentage of dust-only flux at specific wavelengths for each model at 
steUar temperatures T^g > 20000if . Each model is labelled as "t[optical depth] -T [stellar temperature in 1000K]-i[MP/POL]", where MP 
indicates the near mid-plane line of sight (77.5°) and POL indicates the near pole line of sight (12.5°). N/A entries signify cases with 
null denominator fluxes (within the Monte Carlo error). 



The more measurements made at distinct wavelengths the 
better the various physical parameters may be constrained. 

By running MOCASSIN on disc-star systems at vary- 
ing stellar temperatures and optical depths considering ra- 
diative transfer through both the dust and the gas, we have 
shown that irradiation of the gas in discs surrounding high- 
mass stars has a significant effect on the overall shape of the 
SED. These flux differences, particularly in certain wave- 
length ranges, may result in observations of massive YSO- 
disc systems being erroneously fitted if being fitted to "dust- 
only" models, such as the Robitaille et al. (2006) models. We 
cannot quantify this directly, since the models in this paper 
and the Robitaille et al. models make different assumptions 
in the surface density of the disc. In addition, our set of mod- 
els would need to cover a wider range of parameter space in 
order to determine general trends. Nevertheless, we can still 
discuss - qualitatively - issues that are likely to arise if mod- 
eling discs around MYSOs with models that do not take into 
account the effect of the gas: 

• Mid-infrared wavelengths probe the inner region of the 
disc, and therefore are very important for studying disc evo- 
lution, planet formation, etc. Especially for high-mass stars, 
photo-ionizing flux from the central source can ionize and 
evaporate the disc, starting from the inner regions (e.g. Rich- 



ling & Yorke 2000). However, if discs around high-mass stars 
are modeled with radiative transfer models that do not take 
into account the effects of the dust, the inner disc proper- 
ties could be misinterpreted, leading to incorrect conclusions 
about inner disc evolution. 

• The main effect of taking the gas into account on the 
dust in the disc is to lower the incident flux and therefore 
the heating, which causes the far-infrared and sub-millimetre 
fluxes to be reduced. Since this wavelength regime is also 
sensitive to the disc mass, modeling these wavelengths with 
dust-only models would cause the disc mass to be underes- 
timated by more than a factor of two in some cases. 

• If one has a well-sampled SED in the sub-millimeter to 
mm regime, the change in slope due to the free-free emission 
should be unmistakable. However, in most cases, only a few 
fluxes are available at these wavelengths, and it is therefore 
possible that the change in slope would not be seen, but 
instead, the slope of the dust emission could appear flatter. 
This would result in incorrectly modeling the dust proper- 
ties, which would then lead to incorrectly estimating the disc 
mass, disc surface density power, and other disc properties. 

Though the models presented here do indeed indicate 
the need for high-mass YSO-disc systems to be fitted to 
models where gas is considered, there are a number of ele- 
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Figure 2. Ratio of the dust-gas SEDs to the dust-only SEDs for the pole-on viewing angle. 



ments not included in our models that may further affect the 
shape of the SEDs. Firstly, the presence of PAH molecules 
was not considered in the irradiation of these discs. These 
molecules are abundant, ubiquitous and a dominant force in 
the interstellar medium, although little is known about their 
role in the structure and evolution of protoplanetary discs. 
However, the presence of these molecules is inferred through 
strong features in the near to mid-IR. It is possible that the 
inclusion of these molecules in the radiative transfer would 
reduce the difference in flux between the two cases in this 
wavelength range. 



Additionally, the effects of asymmetric accretion onto 
the YSO may result in hot spots on the surface. Without 
further modelling including a detailed accretion spectrum, 
it is difficult to know the effect that these hot spots would 
have on the shape of the SED. In theory however, a gen- 
erally cool star could have regions of its surface with sig- 
nificantly higher surface temperatures, emission from which 
could result in the ionisation of the gas in the disc (perhaps 
in restricted regions) and thus similar changes to the SED 
as those detailed above might be seen. 
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Figure Al. SED plots for mid-plane optical depth Tmp = 0.1. Prom left to right, the SED plots are for the near-pole line of sight 
(12.5°), the near mid-planc line of sight (77.5°), and averaged over viewing angles respectively. From top to bottom, the SEDs are for 
stellar effective temperatures of 5800K to 50,000K. Each panel shoves the SED for the dust-only case (black), and the dust-gas case (red 
to blue). 
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Figure A2. SED plots for mid-plane optical depth Tmp = 1. From left to right, the SED plots are for the near-pole line of sight (12.5°), 
the near mid-plane line of sight (77.5°), and averaged over viewing angles respectively. From top to bottom, the SEDs are for stellar 
effective temperatures of 5800K to 50,000K. Each panel shows the SED for the dust-only case (black), and the dust-gas case (red to 
blue). 
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Figure A3. SED plots for mid-plane optical depth Tmp = 100. Prom left to right, the SED plots are for the near-pole line of sight 
(12.5°), the near mid-plane line of sight (77.5°), and averaged over viewing angles respectively. From top to bottom, the SEDs are for 
stellar effective temperatures of 5800K to 50,000K. Each panel shoves the SED for the dust-only case (black), and the dust-gas case (red 
to blue). 
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